Superfluid stability in BEC-BCS crossover 
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We consider a dilute atomic gas ol two species ol lermions with unequal concentrations under 
a Feshbach resonance. We find that the system can have distinct properties due to the unbound 
fermions. The uniform state is stable only when either (a) beyond a critical coupling strength, where 
it is a gapless superfluid, or (b) when the coupling strength is sufficiently weak, where it is a normal 
Fermi gas mixture. Phase transition(s) must therefore occur when the resonance is crossed. 

PACS numbers: 03.75.Ss, 05.30.Fk, 34.90.+q 



Feshbach resonance |l| has opened up a new play- 
ground for the field of cold trapped atoms. Using this 
resonance, the effective interaction between the atoms 
can be varied over a wide range. In particular, for two 
fermion species with a Feshbach resonance between them, 
the ground state can be tuned from a weak-coupling 
Bardeen Cooper Schrieffer (BCS) superfluid to a strong 
coupling regime where the Fermions pair-up to form 
Bosons which in turn undergo Bose Einstein Condensa- 
tion (BEC) 0,0. 

Though this problem has been under intense theoret- 
ical and experimental p| investigations, almost all 
works thus far are restricted to the case where the con- 
centrations of the two fermionic species are equal. We 
here generalize this study to the case of unequal popu- 
lations of the two species, and investigate in detail the 
thermodynamic stability of this system, in particular the 
question when the uniform state can be stable. 

Studies of fermions with unequal populations or mis- 
matched Fermi surfaces and a pairing interaction have a 
long history. It was studied by Fulde and Ferrell, Larkin 
and Ovchinnikov (FFLO) 0| in the 1960's with relation 
to superconductivity in materials with ferromagnetically 
coupled paramagnetic impurities. It was found that in 
this case the system is likely to have an inhomogeneous 
gapless superconducting phase. Advances in techniques 
of manipulating dilute ultracold atoms have revived in- 
terests in the related problems 0. These studies, in our 
present language, are still restricted to the weak-coupling 
regime. We, however, would extend our analysis to all 
coupling strengths. 

In the "canonical" problem of two species of fermions 
with equal mass (say, spin up and spin down electrons) 
and equal concentrations (thus a single Fermi surface), if 
the cross-species interaction is varied from weak to strong 
coupling, at low temperatures the system would undergo 
a smooth crossover from a superfluid with loosely bound 
Cooper pairs (the "BCS" limit) to one with condensation 
of tightly bound bosonic molecules (the "BEC" limit) 0- 
The situation, however, can be very different if one con- 
siders two species of fermions with unequal concentrations 
(i.e. mismatched Fermi surfaces), even if they have iden- 



tical mass. This can be anticipated because, on the one 
hand, far into the BCS side, the system is basically in the 
FFLO regime Q and therefore must go into a spatially 
inhomogeneous phase. On the other hand, in the far end 
of the BEC side, the system is expected to behave like 
an ordinary (weakly interacting) Bose- Fermi mixture and 
thus has a stable homogeneous phase. Here the bosons 
are the fermion pairs and the fermions are the "leftover" 
unpaired atoms of the majority species. Deep into the 
BEC regime, the size of the Fermion pairs is small and 
the interaction between the bosons and the leftover un- 
paired fermions are expected to be weak. It is therefore a 
very interesting question as to what happens in between. 
This is the question we want to address in this paper. 

For simplicity, we shall assume that the resonance 
is sufficiently wide that the physics reduces effectivel_ 
to a single channel regime. This is probably valid 
for many Feshbach resonances under current experimen- 
tal investigations. Thus, in our calculations, we would 
not invoke explicitly the presence of the "closed chan- 
nel" which leads to this Feshbach resonance. We simply 
model the fermions as interacting through a short-range, 
s-wave effective interaction (dependent on the external 
magnetic field) characterized by the corresponding scat- 
tering length a. 1/a varies from oo for large negative de- 
tuning (closed channel bound state energy much below 
continuum threshold) to — oo for large positive detuning. 

Now we proceed to the details of our calculation and re- 
sults. We consider two fermion species, denoted as "spin" 
t and I, of equal mass m. Because of the unequal concen- 
trations of the two species and the possible existence of 
pairing, it is useful to introduce three fields: the chemical 
potentials [i a (a =f or |) and the pairing field A. We 
shall confine ourselves to zero temperature and general- 
ize the BCS mean field approach of [j . The excitation 
spectrum for each spin is (see e.g. |ll| for details) 
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where £ CT (k) = H 2 k 2 /2m — \i a are the quasi-particle ex- 
citation energies for normal fermions, and — f=L The 
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density of each spin species is then 
d 3 k 



(2tt) 3 



with 
4 



[«k/(30 + vlfi-E-*)} , (2) 

factors 



the coherence 
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function. The equation for the order parameter A reads: 
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We solve equations and self-consistently for fixed 
total density n = 7ij + and density difference = 
n-f — n j,. We shall always take | to be the majority species 
so that 7id > 0. 

It is convenient to introduce the average chemical po- 
tential /j, = + 2 and the difference h = (/x-f — 
Hl)/2 > 0. Then we have 



£ T>i (k) = ^(k)2 + A 2 T h, 



(4) 



where £(k) = ti 2 k 2 /2m — /i. Hence Ei(k) > always. 
From Eq. we get 
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and so the integration is only over the region where 
i?j(k) < 0. In the following, it is useful to note that 
the smallest (or most negative) Ey (k) occurs at £ (k) = 
for (i > 0, where it is A — h, and at k = for \i < 0, 
where it is yj /i 2 + A 2 — h. 

As in the case of equal concentrations, it is convenient 
to express our results in dimcnsionless variables. We 
shall define an inverse length scale fcj? through the to- 
tal density n via kp = {2>ir 2 n) 1 / 3 , and an energy scale 
ep = h 2 k 2 F /2m. We thus write u = y/ e F, h = h/ep, 
A = A/ep, fid = rid/n, and define the dimcnsionless 
coupling constant g = l/^nkpa), which varies from oo 
for large negative detuning to — oo for large positive de- 
tuning. 

We now describe the results of our calculations. We 
first make contact with the BCS-BEC cross-over for equal 
concentrations. The inset of Fig. 1 shows the typical 
behaviors of u, A and h as a function of g for a given 
density difference hd- The behavior of p, or A is similar to 
that in the case of equal concentrations 0. For example, 
fl is large and negative in the BEC limit whereas it is of 
order 1 in the BCS regime. Unlike that case however, 
g has to be larger than a minimum coupling g c in order 
for a finite order parameter A to exist. For g < g c , Eq. 
© requires that A = and the system is in the normal 
state. (For clarity of this inset, we plot only the A ^ 
solutions.) The main Fig. 1 shows h as function of g in 
the intermediate regime (|<?| < 1) for three different fid 
(0.2, 0.5, and 0.8). The horizontal dotted lines indicate 



the normal state in which the gap function A is zero 
(described above). 

The behavior for g < g c is easy to understand. For suf- 
ficently large and negative g, the interaction is too weak 
to produce pairing since the concentrations are unequal. 
The system reduces to a Fermi gas. In this case, the 
chemical potentials are given by \i a — (67r 2 77 (J ) 2 / 3 /(2m) 
which implies /i = [(67r 2 ) 2 / 3 /477i](7i^ 3 + n 2 ^ 3 ) and h = 

[(67r 2 ) 2 / 3 /4m)](n^ /3 -77^ /3 ) (both independent of g). On 
the other hand, for large and positive g (the strong- 
coupling BEC limit), one can show from Eqs. J2|, 
© and JSJ that both h and are large and to lead- 
ing order given by h 2 / (2ma 2 ) . However, (/i + h)/2 = /if 
= (67r 2 7id) 2 ^ 3 /(2TO) <C or h. These expressions sim- 
ply reflect that the system becomes a Bose-Fermi mixture 
with boson concentration nj and free fermion concentra- 
tion rid- 

Notice that the lines (for A ^ 0) of h versus g cross 
each other from small to large hd near g ~ 0.15. For 
g > 0.17, h increases with fid for fixed g. For g < 0.15, h 
decreases as hd increases when the coupling strength is 
fixed. We shall return to these features again below. 

Now we make contact with the superconductivity lit- 
erature. It is helpful here to note that h plays the role of 
an effective external Zeeman field. We plot in Fig. 2 A as 
a function of h for various coupling strengths g. The hor- 
izontal portion of each curve corresponds to rid = 0. In 
this region, h < A and so that (k) > for all k (see Eq. 
(JSJ). The other part of the curve corresponds to rid > 0, 
and exists only in the region h > A (More precisely, for 
larger g where /i becomes negative, this condition should 
read h > \J + A 2 ). For small g (< 0.1), A decreases 
with decreasing h. This solution is the generalization of 
that first discovered by Sarma . We find that this su- 
pcrfluid state corresponds to one where fid increases with 
decreasing h, and hence unstable (to be discussed again 
below). For sufficiently large coupling (g > 0.17), A de- 
creases with increasing h. This state has hd increases 
with h, and satisfies one of the stability conditions to be 
discussed below. 

In Fig. 3 the chemical potential difference h is plotted 
as a function of hd- These results correspond to those of 
Fig. 1 presented in a different manner. Let us explain 
how this graph should be read, with g = —0.1 as an 
example. For hd = 0, h can take any value up to h\ w 0.5 
given by the intersection of the line labelled by g = —0.1 
with the hd — axis. This portion corresponds to the 
line with A being a constant in Fig 2. For < hd < 
0.46 the dependence of h on ti^ is given by the solid line 
labelled by g = —0.1. This line corresponds to the state 
with A ^ but A < h in Fig. 2 discussed above. For 
hd > 0.46, the system enters into the normal state with 
the (/i, hd) relation represented by the dotted lines, given 
by h = 0.5[(1 + h d ) 2 ^ — (1 — h d ) 2 ^ 3 ] (see the discussion 
on Fig. 1 above). Lastly, for hd — 1, h can take any 
value larger than J12 = 0.5 x 2 2 / 3 w 0.79. This is because 
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this line corresponds simply to a Fermi gas with only f 
particles, and h can take any value larger than p so that 
Pi = [p— h)/2 < 0. For g > 0.1, the graph can be read in 
a similar manner except that the dotted line representing 
the normal state is not involved. 

For the uniform superfluid to be stable, two criterions 
must be fulfilled [llL Ha|. First, the susceptibilities ma- 
trix dn a /dp a i can have only positive eigenvalues. One 
can show that this requires that dn^/dh, evaluated at 
constant g, must be positive That is, the plot of h 
versus fid must have positive slope. From Fig. 3, we see 
that for small g (< 0.1), the slope of each curve is always 
negative which indicates the superfluid state is unstable. 
For g greater than about 0.1, the slope of these curves 
change sign at some rid between and 1. In this case 
a stable superfluid state may occur for sufficiently large 
fid- After g > 0.17, the slopes of these curves are positive 
for all rid > and the above stability criterion is satisfied 
for all rid- 

The second stability criterion is that the superfluid 
density p s must be positive [ll|. p s can be evaluated 
as discussed in 0] . The analytic result can be expressed 
as 
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with fci, 2 = [{p T Vh 2 - A 2 )] 1 / 2 and 0{x) is the step 
function. The line p s — is plotted as the dashed lines 
in Fig. 3, with p s < below and p s > above. Thus 
the states correspond to A ^ below this dashed line 
are all unstable. p s < indicates that the system is 
unstable towards a state with spatially varying A and 
hence a state such as FFLO can be more preferable. From 
our results, it turns out that the condition for p s > is 
actually a slightly weaker requirement than the positive 
susceptibility discussed in the last paragraph (though we 
are not aware of any reason why it must be so) . 

We here note that, for fid — > 0, the location where 
p s changes sign is exactly at p = 0. Though this can 
be seen from Eq. @, it is more instructive to return 
to the more basic equation for the superfluid density: 
p s = n + pp. Here p p , the paramagnetic response, is re- 
lated to the backflow of quasiparticles and is given by 
M Pp = -(zkif™dkk*5(E r (k)). For p < the so- 
lution to E-\(k) = exists only when h > \J p? + A 2 
and takes place at k = ki with k 2 /2m = \/h 2 — A 2 + p. 
For rid — > + , ft is just slightly larger than p 2 + A 2 
[see Eq. (0]. k^ is small and hence p v — > + and 
p s ps n > 0. For > 0, E^{k) = happens when 
ft, > A and take place at two k values: k = ki as al- 
ready given in the p < case above and k — k\ with 
k 2 /2m = —\Jh 2 — A 2 + /j. For — > + , ft- is just 
slightly larger than A and the E^(k) = points occur 
near £(fc) w hence dE^(k) / 'dk — > 0. Since and &2 



are finite, p p — > — oo and p s < 0. 

Finally we show our phase diagram in Fig. 4 covering 
the entire BEC to BCS regimes. On the BEC side, with 
g > 0.17, the superfluid state is stable in which both the 
slope of h(fid) and the superfluid density are positive (see 
Fig. 3) . On the upper right of Fig. 4, the pairing order 
parameter A is zero and the system is in the normal 
state. In constructing the boundary of this phase, we 
have simply taken the intersection of the full lines in Fig. 
3 with the dotted lines. A more correct approach would 
involve the Maxwell construction. However, to do this 
we also need to know the equation of state for the non- 
uniform FFLO superfluid state beyond the weak-coupling 
regime. Since this information is not yet available, we 
shall leave this investigation to the future. 

Lastly we discuss the "breached gap" state of Liu and 
Wilczek For this state to exist, one need E^(k) < 
for a region of k that lies between fci < k < where 
k\ > 0. (fci, ki already given in the paragraph before 
last). This is possible only if E^(k) is not monotonic 
with k and hence p > 0. Moreover, since E^(Q) > 0, 

we have \J p 2 + A 2 > h yet h > A. Though there is a 
region of stability with p > (the region near the upper 
right of Fig. 3), we find that it rather corresponds to 
\J p 2 + A 2 < h. (That this is the case near fid ~ 1 is 
obvious, since the line fid — 1 corresponds to h > p and 
A — > 0. Note also on the dash-dotted lines where p = 0, 
we also have \J p 2 + A 2 < h since h > A). Therefore 
gapless excitations exist only near k = k%. Moreover 
this region has n^(k) = 1 for < k < &2 (and n-\(k) = 
ni(k) = v 2 (k) < 1 for k > k^), thus the leftover unpaired 
majority spin-up particles form a rather normal Fermi 
sphere with radius k2, although the energy required to 
create a hole —E^{k) can actually be non-monotonic as a 
function of k, being maximum at an intermediate value 
k = \J1mp (where it is h — A) but not k = (where it is 
h- 



A 2 ). 



In conclusion, we have investigated the stability of a 
fermion mixture with unequal concentrations under a 
Feshbach resonance. We show that, in contrast to the 
case of equal concentrations, there is no smooth BCS- 
BEC crossover. The system is a uniform stable superfluid 
Bose-Fermi mixture only for sufficiently large coupling. 
For weak interactions the normal state is the only stable 
uniform state. The uniform state is unstable for interme- 
diate coupling strengths. Phase transitions must occur 
when the Feshbach resonance is varied between large pos- 
itive detuning and large negative detuning. 
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FIG. 1: Main figure: scaled chemical potential difference 
h/eF = h versus the coupling constant g for three different 
values of n d /n — h d . Inset includes also the results for A and 
p, for rid/n = 0.5. 



FIG. 2: Scaled pair order parameter A/ef = A versus 
h/eF = h for given coupling constants g. 
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FIG. 3: }i/ef = h versus n d /n = h d for constant g's. Full 
lines are for A ^ and the dotted lines are for A = 0. The 
dashed line indicates where p s — 0, with p s > above it. The 
dot-dash line indicates fi = 0, with fi < above this line. 



FIG. 4: Phase diagram for the two Fermion species with 
interspecies coupling g, with the stable uniform phases (white 
regions) shown. No uniform phases are stable in the shaded 
region (except when n d /n — or 1). 



